Jammed frictionless discs: connecting local and global response 
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By calculating the linear response of packings of soft frictionless discs to quasistatic external 
perturbations, we investigate the critical scaling behavior of their elastic properties and non-affine 
deformations as a function of the distance to jamming. Averaged over an ensemble of similar 
packings, these systems are well described by elasticity, while in single packings we determine a 
diverging length scale £* up to which the response of the system is dominated by the local packing 
disorder. This length scale, which we observe directly, diverges as 1/Az, where Az is the difference 
between contact number and its isostatic value, and appears to scale identically to the length scale 
which had been introduced earlier in the interpretation of the spectrum of vibrational modes. It 
governs the crossover from isostatic behavior at the small scale to continuum behavior at the large 
scale; indeed we identify this length scale with the coarse graining length needed to obtain a smooth 
stress field. We characterize the non-affine displacements of the particles using the displacement 
angle distribution, a local measure for the amount of relative sliding, and analyze the connection 
between local relative displacements and the elastic moduli. 

PACS numbers: 45.70.-n, 46.65+g, 83.80.Fg, 05.40.-a 



The jamming transition governs the onset of rigidity 
in disordered media as diverse as foams, colloidal suspen- 
sions, granular media and glasses [3, [|J. While jamming 
in these systems is controlled by a combination of density, 
shear stress and temperature, most progress has been 
made for simple models of frictionless soft spheres that 
interact through purely repulsive contact forces, and that 
arc at zero temperature and zero load [1, d, [f| H, 0, Hi • 
This constitutes possibly the simplest model for jam- 
ming. Moreover, this is a suitable model for static foams 
or emulsions [{J [T(| EH j which also represents a simplified 
version of granular media , ig noring friction (l2l Ha , [l4| 
and non-spherical shapes [l5l [Tfjl . HtI [H| . 

From a theoretical point of view, this model is ideal 
for two reasons. First, it exhibits a well defined jamming 
point, "point J", at confining pressure p = 0, and in the 
limit of large system sizes, the jamming transition oc- 
curs for a well-defined density <j) — <p c , which has been 
identified with the random close packing density [|[ . At 
this jamming point, the system is a disordered packing of 
frictionless undeformed spheres, which is marginally sta- 
ble and isostatic, i.e. its contact number (average number 
of contacts per particle) z equals Zj SO = 2d in d dimen- 
sions. Second, in recent years it has been uncovered that 
the mechanical and geometric properties of such jammed 
packings of soft spheres close to point J exhibit a num- 
ber of non-trivial power law scalings as a function of 
A(f> := cj) — <j) c . These scalings illustrate the unique char- 
acter of the jamming transition [1, 0, H, @, 0, Hi ■ 

Approaching point J from the jammed side, the most 
important scaling relations are as follows. (1) The excess 
contact number Az := z — z iso scales as (A0) 1 / 2 [H, 0, 
1@. (2) The ratio of shear (G) and bulk (K) elastic 
moduli vanishes at point J as G/K ~ Az (3) The 
density of vibrational modes exhibits a crossover from 



continuum like behavior to a plateau at a characteristic 
frequency lo* , which vanishes at the jamming point: lu* ~ 
Az [1, [a • O ne can associate a diverging length-scale 
£* with this crossover, which then diverges at point J: 
I* ~ 1/Az @. 

Recently, we uncovered an additional non-trivial scal- 
ing near point J when identifying the degree of non- 
affinity Decomposing, for linear deformations of 
jammed systems, the relative displacement uy of neigh- 
boring particles i and j in components parallel (u||) and 
perpendicular (u±) to ry, where ry connects the centers 
of particles i and j, we find that the ratio u±/u\\ diverges 
near point J. More precisely, defining local displace- 
ment angles ay via tan ay = ux,ij /uu^j , we find that 
the displacement angle distribution P(a) peaks around 
a = 7r/2, with the peak height diverging near point J. 
Hence, close to point J, the non-affinity is such that par- 
ticles predominantly slide past each other. 

The nature of the non-affinity of particle displacements 
ties in with the question to what extent these disor- 
dered solids close to the jamming transition can be de- 
scribed using continuum elasticity. Recent work on gran- 
ular materials has shown that this is possible, using the 
proper definitions of stress and strain tensor (f9j . and 
provided one probes the system at a large enough length 
scale 0, [2l|, . However, it has remained unclear if this 
length scale is to be identified with the crossover length 
scale £* that was introduced in the interpretation of the 
density of vibrational states @. In addition, despite the 
large body of numerical and experimental work on the 
problem [2I, [H, [H, [H, S3, H, , a systematic analy- 
sis of the relations between the bulk clastic moduli, the 
stress fields in response to local perturbations, and the 
distance to the jamming transition appears to be lacking. 

In this paper, we describe the applicability of elastic- 
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ity theory for jammed packings, and elaborate on our 
earlier work on the displacement angle distributions @. 
Although all our numerical studies are restricted to fric- 
tionless spheres in two dimensions, we expect essentially 
all concepts and conclusions to carry over to three di- 
mensions as well. 

In section I we detail our linear response methodology, 
and introduce our notation. In section II we focus on the 
response of jammed systems to local and global forcing, 
and show that this response, suitably coarse grained and 
averaged, can be described by linear elasticity. Earlier 
direct numerical simulations [j, |30| have established the 
scalings of the shear modulus G and bulk modulus K 
with distance to jamming, and in particular have shown 
that their ratio G/K vanishes at point J. Our linear re- 
sponse calculations reproduce these findings. 

In section III we address the issue of the length scale 
beyond which elasticity theory can be applied to the sys- 
tem. First, we determine this length directly from the 
response to a local perturbation, as the scale up to which 
the response is dominated by spatial fluctuations, and 
identify it with I* ~ 1/Az @. In addition, we show that 
the coarse graining length, needed to obtain a smooth 
response in a globally deformed system, is proportional 
to the same length £*. 

In section IV we characterize the non-affine nature of 
the response to bulk deformations, and elaborate on the 
scaling of the displacement angle distribution P(a). We 
discuss the connection between this scaling and the form 
of the expression for the changes in clastic energy in lin- 
ear order: the opposite sign of the contribution from the 
parallel (u\\) and perpendicular (it j_) components of the 
local displacements naturally leads to a delicate balanc- 
ing act for the case of an overall shear deformation, while 
compression leads to a more convoluted picture [3l| . 



I. LINEAR RESPONSE 

All numerical results presented in this paper concern 
the quasistatic linear response of a certain starting config- 
uration to a small perturbation. The approach is there- 
fore explicitly a two-step one. First, we prepare a two- 
dimensional packing of frictionless polydisperse discs at a 
given density or confining pressure, using a Molecular Dy- 
namics simulation. The resulting packings are then ana- 
lyzed by studying their response to small perturbations. 
Response studies have been done previously by using MD 
also for the perturbed system [30], in which case the full 
dynamics are taken into account. Another method that 
has been used is the quasistatic approach of minimizing 
the potential energy, while ignoring inertia 0]. We here 
take the even simpler approach of explicitly focussing on 
the linear response equations. This requires solving just 
a single matrix equation for each numerical experiment, 
which makes it numerically cheap. 

In the following subsections we describe the procedure 
to generate our packings [HI, [32|, and recapitulate the 



derivation of the linear response equation (8l. l33l. [33 from 
an expansion of the elastic energy of the system 0, [H[ , 
both for completeness and to introduce the necessary no- 
tation. 



A. Packing Generation 

We prepare our packings using a molecular dynamics 
simulation with round discs in two dimensions. The discs 
interact via the three-dimensional Hertzian potential, 

i,, r t) 5 nj - dij , a) 

{ nj > dij 

where i, j label the particles, dij is the sum of the particle 
radii Ri + Rj , and r»j is the distance between the particle 
centers. The energy scale depends on the radii and ef- 
fective Young moduli of the particles, sec Appendix I A II 
The quantity between brackets is the dimensionlcss over- 
lap, 6ij = 1 - rij/dij. 

Using the 3D Hertzian potential makes the system a 
quasi-2D packing of discs with round edges. We use zero 
gravity to have a homogeneous pressure, and the radii 
of the discs are drawn from a flat distribution between 
0.4 < R < 0.6, thus creating a polydispcrsity of ±20% 
around the average particle diameter d = 1 (our unit of 
length) to avoid crystallization. The simulation starts 
from a loose granular gas in a square box with periodic 
boundaries, which is compressed to a target pressure p 
by changing the radii of the particles while they are mov- 
ing around. The radii are multiplied by a common scale 
factor r s , which evolves in time via the damped equation 
r s = -4w f s - Lul[p(t,r s )/p - l]r s , where uj ~ 6 • 10~ 2 , 
p(t, r s ) is the instantaneous value of the pressure and p 
the target pressure. This ensures a very gentle equili- 
bration of the packings. Energy is dissipated through 
inelastic collisions and a drag force which slows down 
the particles. The simulation stops when the accelera- 
tions of all grains have dropped below a threshold which 
is 10~ 6 (/) in our reduced units. This way we generate 
packings in mechanical equilibrium for pressures ranging 
from p = 10~ 6 to p = 3 ■ 10 -2 , in units of the modified 
Young modulus of the individual grains (3(|. This corre- 
sponds to a range of contact numbers from z = 4.05 to 
z = 5.87. For one particular type of calculation we use 
packings which are only periodic in the x-direction, and 
have hard walls on top and bottom. These are gener- 
ated by compressing the system using a hard piston, as 
described in Ref. [321 ]. 

Because there is no gravity in our simulations, at the 
end of the simulation there will usually be particles with- 
out neighbors, or due to numerical precision effects, par- 
ticles with fewer neighbors than is needed to make them 
rigidly connected to the rest of the packing. These rat- 
tlers or floaters do not contribute to the rigid backbone 
of the packing. They are removed from the system when 
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FIG. 1: The distance to the jamming transition is set in our 
packing preparation by the pressure p. Other parametriza- 
tions are (a) the distance to isostaticity Az = z — z\ ao w 
9-p ' 37 and (b) the typical dimensionless interparticle overlap 
5 w 1.4 -p . The solid lines represent these empirical power 
law fits. Diamonds represent periodic packings, and squares 
are for packings with hard walls on top and bottom. 



determining the contact number z or calculating its linear 
response to small perturbations. 

To express the distance to the jamming point, we will 
use the pressure p and typical relative overlap 5 inter- 
changeably; the pressure (rather than the density) is 
most conveniently set in the numerical procedure to gen- 
erate the packings, while the overlap enters in the scaling 
relations at the particle scale that will be discussed in sec- 
tion [TV] Wc do not use Acj) = (j> — <f> c , because in finite 
systems (f> c is a number that would have to be determined 
for each packing separately. However, close to point J wc 
have determined that A(f> ~ S. The conversion between 
the different parametrizations is given in Fig. [TJ 



B. Linear response equation 

We calculate the response of a packing to a mechani- 
cal perturbation, which can be either an external force or 
a deformation of the periodic box. The response in gen- 
eral involves displacements of all particles in the packing, 
which wc analyze in the linear regime. 

The total potential energy of the system is a sum over 
all pairs of interacting particles, 



(i,3) 



(2) 



where we assume we are dealing with a central poten- 
tial, only depending on the distance between the parti- 



cles 



= r., - r, ; 



The change in due to 



displacements of the particles is 



A? 



(3) 



where mm is the relative displacement along the line con- 
necting the centers of the particles and u± is the rela- 
tive displacement perpendicular to that. In the linear 
response regime the change in energy is expanded to sec- 
ond order in tin and uj_ as 



AE 



*5> 



fk 



(4) 



Here — — dT4,/dr„ and fey = d 2 Vy/dr|j. For all 
(power law) interactions that are reasonable models for 
foams or granular media (linear repulsion, Hertzian re- 
pulsion), both the initial force / and the stiffness k are 
nonnegative (see below). There are no terms linear in 
the displacements because the starting configuration is in 
mechanical equilibrium, which makes these terms sum to 
zero. The w|-term represents the change in bond length. 

The w^-term is only there if there is a pre-stress or ini- 
tial force and captures the lowering of the energy due to 
transverse displacements of the particles [13] ■ 

The interaction potential for the particles that make up 
the packings used in this paper is a finite-range, purely 
repulsive 1 potential of the form \ ' 



8fj, where 



(5) 



is the dimensionless virtual overlap of the particles and 
Ri and Rj are the radii of particles i and j. The exponent 
a has the value a = 5/2 in our packings, representing the 
three-dimensional Hertzian interaction. For any poten- 
tial of power law form, the prefactor of the second term 
in equation ^ can be written as &ij/{a — 1). The closer 
to point J, the smaller these dimensionless overlaps, so 
this prefactor will be small close to point J. However, this 
does not allow us to ignore this term, because, as we will 
see, the typical u± arc going to be much larger than the 
typical it || , in the limit of approaching the jamming tran- 
sition: the two terms in Eq. ([4]) become, in some cases, 
of the same order. 

Writing the change in energy in the independent vari- 
ables of the problem, the displacements of the parti- 
cles, defines the dynamical matrix M: 



AE 



(G) 



where Ai is a dN x dN matrix (d being the spatial di- 
mension and N the number of particles). The indices a, (3 
label the coordinate axes and we use the summation con- 
vention. The dynamical matrix contains all information 
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on the elastic properties of the system. It can be diago- 
nalized to study the vibrational properties [1, [TH, [35[ , but 
it can also be used to study the response of the system 
to an external force (defined on each particle) [1, [H, [H| : 

M ijta0 u h p = ffg . (7) 

The dynamical matrix M. is very sparse for large systems, 
because the only nonzero elements are those for which i 
and j are in contact and those for which i — j. There- 
fore we can compute the response efficiently by using the 
Conjugate Gradient Method This is essentially an 

iterative procedure that minimizes WM-ij^pUj^ — ff^\\- 
O'Hern et al. [1] also studied the quasistatic response of 
granular systems to a global shear or compression Q us- 
ing a Conjugate Gradient Method, but in their case the 
quantity to be minimized was the potential energy. The 
difference is therefore that we use Conjugate Gradient 
only as an efficient method to study small deviations from 
a stable starting configuration in linear reponse, while 
O'Hern et al. used it for various situations where they 
needed to find the nearest potential energy minimum. 

C. Forces and stresses 

Solving the linear response equation ([7]) yields the dis- 
placements Ui of all particles. From this we calculate 
the local relative displacements u±^j and wiuj, and the 
change in contact force A/ij using 

= cos(/)ijUij tX + sin <f>ijUij !V , (8) 
u±,ij = cos^y Uy^ - sin^y Uy^ , (9) 
A/ij = -ku\\ ti j , (10) 

where fyj is the angle which the bond vector makes 
with the x-axis. 

To get from discrete forces to a continuum stress field, 
we apply the local coarse graining procedure developed 
by Goldhirsch and Goldenberg fig ]: 

Acr Q/3 (r) = ^2 A/y, a '"y, / 3 / ds $ (r - v + sr y ) . 

(11) 

It has been shown that this procedure gives local stress 
fields that do not depend strongly on the coarse grain- 
ing length chosen — consistent results were obtained for 
length scales down to a single grain diameter [22|. In 
this paper, we use a Gaussian coarse graining function $ 
with a width £cg equal to the average particle diameter, 
£cg = 1 in our units, i.e. 

$(r) = — e -l r l 2 /2. 

If, instead, one uses $ = l/V, Eq. (TTTj) reduces to 
the well-known expression for the average of the contact 
stress tensor over the entire system. 



The precise form of Eq. (jlip is not crucial for the work 
presented in this paper, because we will only use this 
coarse graining procedure for the stresses, not for the 
strains. It is the requirement of emerging linear elasticity 
from coarse graining both stress and strain that led to 
this particular form in Ref. 19]. It should also be noted 
that, in principle, the expression for Act contains both 
terms ~ A/ r and ~ / Ar — we ignore the latter since 
these terms are negligibly small close to the jamming 
transition. 



II. MACROSCOPICS AND CONTINUUM 
ELASTICITY 

There are three questions that will be discussed here 
concerning the macroscopic aspects of the linear response 
of granular packings. First of all, under what conditions 
can the system's response to external forcing be described 
using continuum elasticity? Secondly, do we recover, in 
the linear regime, the scaling of the elastic moduli with 
contact number that was obtained in direct numerical 
simulations 0? Finally, is there a difference in the values 
of the elastic moduli as calculated from the response to 
global forcing, using bulk compression or shear, and local 
forcing, applying an external force on a single particle? 

The work by Goldenberg and Goldhirsch 0, H], HH, 
[22| has made clear that stress and strain tensors are well- 
defined down to even less than the scale of a single grain, 
using the coarse graining expressed by Eq. (jlip . For 
large enough systems, they find elasticity to provide a 
good description of the response of granular media. In 
this section, we analyze the response of granular packings 
as a junction of distance to the jamming transition, and 
confirm that the linear response of granular media, when 
averaged over an ensemble of similarly prepared pack- 
ings, is well described by continuum elasticity. To do so, 
we study the response to both global and local perturba- 
tions, as described in sections III Al and III B[ respectively. 
A fit of the ensemble averaged stresses and displacements 
provides the elastic moduli of the packings, which we find 
to scale consistently with earlier results |3| . 



A. Bulk deformations 

The conventional way to extract elastic constants of a 
packing is to apply a compression or shear deformation to 
the entire system. In packings with periodic boundaries 
this is done by imposing a relative displacement on all 
bonds that cross the boundary of the periodic box, fol- 
lowing a procedure that is similar to the Lees-Edwards 
boundary conditions employed to simulate uniform shear 
flows |39j. For example, to impose a globally uniform 
compressional strain e xx = e yy = e, all terms in the en- 
ergy expansion of equation ^ that represent a bond that 
crosses the periodic x-boundary are changed such that 
each occurrence of Uj — is replaced by Uj — U; ± eL x x, 



FIG. 2: (color) Force response networks for a point loading with pressure as indicated. Blue (red) lines indicate an increase 
(decrease) of the contact force, the thickness indicating the amount (the thickness of the black border around each panel 
corresponds to 1/15 of the loading force on the center particle). The particles themselves are not drawn. The figures show only 
the part of the packing close to the forced grain, and contain about 3500 of the 10 000 particles. 



where the ±-sign is given by the sign of Tij iX . The y- 
boundary is treated analogously. Shear deformations are 
applied in the form of a pure shear, i.e., by having a 
displacement in the y-direction imposed on all bonds 
that cross the x-boundary and a displacement in the in- 
direction on all bonds that cross the y-boundary. 

When we write the linear response equation (J7]) from 



the energy expression for the deformed system, the terms 
proportional to e end up on the right hand side of the 
equation and act like an effective / ext . The response of 
the system to this shape or volume change of the periodic 
box is again calculated by solving equation ([7]) for this 
effective external force. The elastic moduli then follow 
from the resulting change in stress tensor according to 
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equation (jlip with the trivial coarse graining function 
$ = l/V. The bulk modulus is extracted from a uniform 
comprcssional strain e xx = e yy = e through 



K 



^yy 
a,. 



2,c aa 



4e 



(12) 



and the shear modulus from a uniform shear strain e xy 
e through 



G 



>xy 



2e 



u xy 

~2e~ 



(13) 



The results will presented and discussed in section III CI 
together with those obtained from the local point re- 
sponse. 

B. Point response 

Now let us determine the linear response of packings 
of ./V = 10 000 particles to a force in the y-direction, ap- 
plied to a single particle. The packings are periodic in 
the x-dircction and have hard walls on top and bottom 
to carry the load. The confining pressure p used to create 
the packings ranges from p = 10~ 6 to p = 10~ 2 , corre- 
sponding to contact numbers ranging from z = 4.05 to 
z = 5.41. We have 10 different packings at each pressure, 
and use each one about 20 times by applying a point force 
to different particles, all of which are closer than 0.1 par- 
ticle diameter to the horizontal line through the center of 
the system. Examples of the resulting response networks, 
depicting the changes in the contact forces, are shown in 
Fig. [21 These pictures immediately illustrate to the eye, 
that as the jamming point is approached, both the range 
and the magnitude of the fluctuations increase. Quanti- 
fying this behavior is the goal of section IIIII — here we 
focus first on the ensemble averaged response. 

In order to allow comparison to continuum solutions, 
we first calculate for each force response network the asso- 
ciated stress response fields, by the coarse graining pro- 
cedure in Eq. CP. See Ref. [3 for details. We then 
calculate the ensemble average of these stress response 
fields, and results of this procedure are shown in Fig. [3] 
The continuum solution that we are comparing the gran- 
ular point response to, is obtained from the static Navier- 
Cauchy equation (40| , which is given by 



GAu + A'V(V • u) 



(14) 



in a two-dimensional formulation of elasticity theory. 
Here K and G denote bulk and shear modulus. This 
equation is the direct equivalent of the linear response 
equation ([7]), because solving it yields displacements in 
terms of external forces. We solve equation (fT4"|) with 
pxt e q ua i t a un it point force in the y-direction; see 
Appendix [Bl for details. The resulting stress field <r Q /3(r) 
only depends on the ratio K/G of elastic constants, since 
the overall scale drops out when relating an imposed force 
to a resulting stress. Therefore, there is only one free 




FIG. 3: Stress response fields for the linear response to a point 
force. The greyscale images represent the ensemble averaged 
granular stress response field; the thick contours denote the 
fitted continuum stress field. For a yy thin contours for the 
granular stress response field are included as well. The com- 
ponents of the stress tensor are plotted for (a) p = 10 -2 , (b) 
p= HT 4 , (c) p= 1(T 6 . 



parameter when fitting the stress field of the granular 
system, equation (jlip . to the continuum expression. As 
described in Appendix [Bl we use the the effective Poisson 
ratio v = (K — G)/ (K + G) . To determine both mod- 
uli we need a second fitting parameter, using a fit to the 
displacement field. In particular, we fit the the average 
y-displaccmcnt of all particles in a strip at height y, 



H grll n(y) := {uLy) yiKy , 
to its continuum counterpart: 



H(y) := ^ 



W2 



irr/2 



u y (x,y)dx 



L v -2\y\ 
AL X (K + G)' 



(15) 



(16) 



which is obtained by taking f oxt = S(x)S(y)y and inte- 
grating equation (|14[) over x, leaving a standard Green's 
function problem with the boundary condition that u y 
vanishes at the top and bottom wall, i.e., at y = ±L y /2. 

Figure[3]displays the ensemble averaged stress response 
fields Aa a f3(r) (equation (fTTj) ') and the fitted continuum 
stress fields for p = 10~ 2 , p = 10~ 4 , and p = 10~ 6 . The 
grey scale and the contour values are chosen for each 
tensor component separately, but for each component are 
the same for the different pressures. The fit is very good, 
especially considering the fact that a tensor field with 
three components is fitted with only 1 parameter. For the 
yy-component we also include contours for the numerical 
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FIG. 4: The fitting procedure used to obtain K + G: The 
average vertical displacement of particles at height y, equa- 
tion (|15[) , is fitted to the triangle-shaped function H(y), given 
by equation (|16p . The functions are rescaled by the fitted 
K + G and vertically offset for clarity (note that H{y) = for 

V/Ly = ±\). 



data, to give a more quantitative view of the fit. For the 
xx- and ^-components the numerical data is too noisy 
for this to be useful, especially at the lower pressures. 

The fits of the average vertical displacements H gT i in (y) 
are shown in Fig. [4j Again a good fit is obtained with 
only one parameter, with only a little bit of noise for the 
lowest pressure packings. Combining the results of the 
two fitting procedures yields the elastic moduli — the 
results are presented in the next subsection. 



C. Results 

The clastic moduli resulting from the methods ex- 
plained in the previous two subsections are collected in 
Fig. [5h.. The squares and diamonds represent the bulk 
(K) and shear (G) modulus, respectively, calculated from 
the linear response to a bulk compression or shear. These 
have been obtained via ensemble averages over 100 pack- 
ings of N = 1000 particles, for pressures ranging from 
p = 5 ■ 10~ 6 to p = 3 • 10~ 2 (in units of the modified 
Young modulus of the grains). This pressure range cor- 
responds to a range in contact number from z = 4.10 to 
z = 5.87. The cross and plus signs indicate the moduli re- 
sulting from fitting the stress response to a point force to 
continuum elasticity. As mentioned in the previous sub- 
section, the packings used for the point response calcula- 
tions were prepared at pressures ranging from p = 10~ 6 
to p = 10~ 2 , corresponding to contact numbers from 
z = 4.05 to z = 5.41. 

From Fig.[3]one sees that the data is well described [54| 
by scaling relations of the form 



K 
G 



0.36±0.03 



0.70±0.08 



(17) 
(18) 



These are consistent with the scalings K ~ P 1 ^ 3 an d 
G ~ p 2 / 3 which on the basis of previous work Q are ex- 



(a) 10° 
10 1 

. - lO" 2 

lO" 3 

10" 4 
10" 5 



K 



10- 



G 



P 



io- 



(b) 



1 



o.i 



X 



.0 



X 



)( 



0.01 



0.01 



0.1 



10 



Az 



FIG. 5: (a) Bulk modulus K and shear modulus G as a func- 
tion of pressure. The squares (K) and diamonds (G) are 
obtained using the bulk response described in section III Al 
The crosses (K) and plus signs (G) follow from the fits of the 
point response stresses and displacements discussed in sec- 
tion HlBl The error bars on the bulk response data indicate 
the intervals spanned by the median 50% of the data; the 
actual standard error of the mean is much smaller than the 
symbol size. The error bars on the point response data are 
error estimates from the fitting procedure. We attribute the 
fact that the point response result for G at the smallest pres- 
sure deviates from the scaling behavior to the fact that our 
fitting to elastic continuum behavior is very insensitive to the 
value of G in the limit G <SC K. (b) The ratio of elastic moduli 
G/K scales approximately as Az. 



pected for our 3D Hertizan contacts. The above scaling 
behavior already shows that we can consistently describe 
the response in terms of continuum elasticity. 

Alternatively, we can plot the ratio of the elastic mod- 
uli as function of Az, and find that G/K scales as Az 
(Fig. . This result is expected to be independent of 
the force law — in all disordered packings that have been 
checked, the bulk modulus K is proportional to the con- 
tact stiffness k, while the shear modulus G behaves pro- 
portional to kAz [1, H3, 5H. The scaling of the shear 
modulus has been described as anomalous 0, [ID, [30| . 
because earlier effective medium theories had predicted 



FIG. 6: Comparison of the stress response to a point force 
of a packing close to the jamming transition (p — 10~ 6 ), (a) 
calculated with the energy expansion as stated in Eq. Q, 
and (b) calculated without the u^-term, which corresponds 
to ignoring the forces present in the system before the point 
force was applied. Thick curves are the continuum fit, thin 
curves the numerical contours, as in Fig. [3] 

K ~ G ~ k [l2|, IHj]. However, from the perspective of 
rigidity percolation on a random network, for which all 
clastic moduli are proportional to Az, the compression 
modulus can be described as anomalous. We have re- 
cently explored this issue extensively in the context of 
harmonic networks [3lj ] - 

To our knowledge, our calculations are the first to show 
that the average large scale response to both a local force 
and a global deformation or force is consistent with con- 
tinuum elasticity theory with the same clastic moduli. 
We thus conclude that the linear response to a local per- 
turbation, averaged over an ensemble of packings, is well 
described by linear elasticity. To what extent and on 
what length scale the response of a single packing cor- 
responds to elastic behavior will be the subject of sec- 
tion nj 

Let us close this section by pointing out the effect of 
the second term in the energy expansion, Eq. which 
is proportional to Uj_. This term contains the influence 
of the forces that the contacts carry before applying the 
external point force, and is therefore referred to as the 
pre-stress term (37|- This term is strictly negative in a 
system of purely repulsive particles, and therefore tends 
to destabilize the system, to enhance the perpendicular 
(sliding) motions of particles, and to lower the energy as- 
sociated with deformations. To study its effect, we also 
performed a calculation in which we left it out. The 
destabilization can be seen from comparing Fig. [5^ with 
Fig. [lb: the spatial fluctuations from the continuum elas- 
ticity stress fields are much smaller if we leave out the 
pre-stress term. The fact that the elastic energies are 
lower when the term oc u\ is included in the energy ex- 
pansion can also be shown by calculating the elastic mod- 
uli with and without this pre-stress term. The resulting 
shear modulus without this term is higher than when 
this term is included. We will come back to the relative 
importance of the two terms in the energy expansion in 



FIG. 7: (color) Force response networks for inflation of a sin- 
gle particle, with p as indicated. Blue (red) lines indicate 
positive (negative) changes in contact force, the thickness in- 
dicating the amount, which is drawn for a 2% increase of 
the particle diameter. The thickness of the border around 
the panel corresponds to a change in force of 1/4500 (left), 
and 1/60000 (right), respectively, which is needed because the 
shear modulus is much lower at the lower pressure. The par- 
ticles themselves are not drawn. The figures show only the 
part of the packing close to the inflated grain, and contain 
about 1800 of the 10 000 particles. 

section IIVI 

III. CRITICAL LENGTH SCALE 

In the previous section we have seen that the ensem- 
ble averaged response of a jammed granular medium can 
be described using continuum elasticity. On the other 
hand, the disordered nature of the packing has a strong 
effect on the displacements and forces in individual real- 
izations, especially in systems close to the jamming tran- 
sition. The question is whether we can think about the 
role of disorder as relevant on small length scales, but 
sufficiently smeared out on long length scales. We will 
probe this question by locally forcing the system, and 
studying the fluctuations away from the ensemble av- 
eraged response as a function of length scale. We will 
find a length scale I* which indicates above what length 
one can consider the system a continuum medium. This 
length scale is proportional to 1/Az, which implies that 
it diverges as we approach the jamming transition. 

The length scale £* was introduced earlier by Wyart to 
describe the excess of low frequency modes in the den- 
sity of vibrational states of these systems @, where it 
can be derived as follows. If a system of dimension d is 
to be described as a continuum medium, it should keep 
its properties when cut open and split in two parts. In 
particular, if we cut out a circular blob of radius £ of a 
rigid material, it should remain rigid. The rigidity (given 
by the shear modulus) of jammed granular materials is 
proportional to Az, the density of excess contacts over 
the minimum number of contacts required to be mechan- 
ically stable. The circular blob has of the order £ d Az of 
such contacts. To cut it out, however, we had to break 
the contacts at the perimeter, of which there are of order 
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zl x . If the number of broken contacts at the edge is 
larger than the number of excess contacts in the bulk, the 
resulting blob is not rigid, but floppy (see Appendix [Cj) . 
The smallest blob one can cut out without it being floppy 
is obtained when these numbers are equal, which implies 
that it has radius £* ~ z/Az. Close to the jamming 
transition, z is essentially constant, and so one obtains 
as scaling relation that [3] 



(a) 



I 

~Az 



(19) 



So far this length scale £* has not been observed di- 
rectly. What has been observed is a crossover frequency 
lu* in the density of vibrational states, marking the 
lower end of a plateau of excess states. The vibrational 
states at uj < co* have been interpreted as ordinary plane 
wave-like modes d, @. Using the speed of sound one 
can therefore translate the crossover frequency into a 
wavelength, which scales as At ~ 1/y/Az for transverse 
(shear) waves and as Al ~ 1/Az for longitudinal (com- 
pressional) waves. At has been observed in the spatial 
structure of the vibrational modes by Silbert et al. [J], 
but it is Al that coincides with the length scale £* derived 
above. Note that the derivation of I* involved neither 
shear nor compression waves and therefore the corre- 
spondence of £* and Al is not obvious a priori. Below we 
present the first real-space observation of £* , in the fluc- 
tuations of the force response to a localized perturbation. 
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FIG. 8: Identification of the length scale I* . (a) The fluctua- 
tions in the contact forces, measured by h(r), which is defined 
in equation (|20[) . are larger in the response of packings at low 
pressure, (b) The same data on a double logarithmic plot, 
showing that the tail decays as a power law, approximately 
1/r 1 ' 6 . (c) The data collapses when we rescale the r-axis by 
Az, signaling a length scale that scales as £* ~ l/Az. (d) The 
relative fluctuations (fluctuations divided by the average value 
of the contact forces). 



A. Observation of £* in inflation response 

The signature of the length scale £* can readily be ob- 
served in the point force response networks (see Fig. [2]): 
The lower the pressure, and hence the smaller Az, the 
larger the scale up to which the response looks disordered. 
To study this effect quantitatively, we calculate the re- 
sponse to an inflation of a single central particle (Fig. [7J 
instead of directional point forcing (Fig. [2]). This allows 
us to probe a natural length scale of the system, as we 
expect a crossover between the behavior for small and 
large r: far away from the inflated particle we expect a 
smooth response with radial symmetry, for which Aa rr 

A G 

Aa rr w — , 

while nearby, the disorder dominates the response. As 
we will show now, the crossover length can be identified 
with £*. 

Examples of the changes in contact forces Af in re- 
sponse to the inflation of a single particle are shown in 
Fig. for pressures p = 4 • 10 -3 and p = 5 • 10~ 5 . For a 
fixed change in radius of the central particle, the scale of 
Af is strongly influenced by the stiffness of the few con- 
tacts of this particle, leading to large fluctuations in the 
amplitude of Af. Since we are mainly interested in the 



spatial structure rather than the values of A/, we first 
normalize the force response Af. To do so, we fit the 
change of each local radial stress Aa rr to the continuum 
field (with a correction for the periodic boundaries), 



and define the normalized response A/ n as Af/Gat- In 
exceptional cases, the fitting parameter Ggt is anoma- 
lously small. Inspection of the response networks in 
which this happens reveals relatively large "soft spots" , 
regions where the rearrangements of the particles are 
large. We leave an analysis of these soft spots and their 
relavance to the future; for the present analysis we limit 
ourselves to noting that these soft spots sometimes lead 
to very bad fits of the stress field, and discard the 1% 
worst- fitting response networks. 

In view of the radial symmetry, on average, we study 
the fluctuations in the response at a distance r from the 
inflated grain. More precisely, we calculate the root mean 
square fluctuations of the radial component of the nor- 
malized force response, A/", through all contacts at that 
particular distance r: 



h(r) = ([Af?(r) - (Af ? (r))] 2 ) . (20) 

Here the average (•) is taken over all bonds that intersect 
the circle of radius r centered around the inflated particle. 
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Note that h[r) is not a simple correlation function in the 
ordinary sense: it simultaneously involves all contacts at 
a distance r from the center and cannot be expressed in 
terms of a single n-point correlation function. We cal- 
culate this function for packings of N = 10 000 particles 
and average over 356 response networks at each value of 
p. The resulting h(r) are shown in Fig. [5^- As expected, 
the fluctuations arc larger and longer ranged for packings 
at smaller pressure. Figure [8Jd shows the decay on a dou- 
ble logarithmic plot: it appears that for large r the tail of 
h(r) goes roughly as 1/r 16 while the small r behavior has 
a smaller slope that varies with p. The crossover scale 
between these regimes is proportional to the length scale 
£*, as is illustrated by the data collapse that is obtained 
when h(r) is plotted as function of rAz in Fig. Eh. One 
might expect the fluctuations to decay as r~ 2 , following 
the decay of the stress. We attribute the difference in 
exponent to finite size effects. 

In Fig. [SJi we plot the relative fluctuations with re- 
spect to the average radially transmitted force (Af^(r)) 
(Fig. [5]i). These are nearly constant (<~ r 4 ) after 
rAz f=a 6, which can serve as a choice of prefactor for 
the crossover length: 




FIG. 9: (color) Force response networks for a global shear 
deformation of packings of 1000 particles at high (left) and 
low (right) pressure. The thickness of the border around the 
panel corresponds to a change in force of 0.5 (left), and 80 
(right) per unit strain, respectively, which is needed because 
the shear modulus is much lower at the lower pressure. Note 
that in the high pressure case nearly all bonds that have an 
angle between and n/2 with the z-axis are red (decreased 
force) and nearly all bonds in the other diagonal direction are 
blue (increased force), consistent with the observation that 
the displacement fields for this shear response are very close 
to affine (Fig. 1110 . The response network of the low pressure 
packing is much more disordered. 



For large r the relative fluctuations do not go to zero, 
which indicates the system is not self-averaging. The 
length scale I* thus marks the distance at which the 
relative error stops growing. Note again that the r in 
this analysis refers to the distance from the perturbation 
where the flucuations are measured, and that the fluctu- 
ation analysis does not involve a coarse graining length. 
The fluctutations, at any r, are measured on the scale 
of single contacts. This contrasts with the next subsec- 
tion, where we will relate I* to the coarse graining length 
needed to obtain a smooth response. 




0.1 1.0 10.0 0.1 1.0 10.0 



FIG. 10: (a) Inhomogeneity H of the response as a function of 
coarse graining length £cg, for various pressures. The lines in- 
dicate the l/£cG-behavior and the error bars denote the 20% 
median values of our dataset of 90 packings at each pressure, 
(b) The same data, with the horizontal axis rescaled by Az 
to collapse the data. The line again denotes H ~ (£cgAz) - . 



B. Inhomogeneity of global response 

The results for the crossover length suggest that to ob- 
tain a smooth stress response field by coarse graining A/, 
one would need to use a coarse graining length propor- 
tional to £*. To test this explicitly, we now study the 
stress fluctuations in packings under an overall shear as 
a function of the coarse graining length £cg ■ 

The starting point are force response networks for 
packings of 10 000 particles that we obtain in linear re- 
sponse to shear (Fig. [5]). The behavior is consistent with 
what we observed for point forcing: the characteristic 
length scale of these force fluctuations appears to grow 
when approaching the jamming point. To quantify this, 
we calculate the change in shear stress Aa xy in 64 points 
using equation (|11[) with a Gaussian coarse graining func- 
tion. The (relative) standard deviation of the resulting 64 



values is a measure of the inhomogeneity of the response: 



^ CG ^(A^)v( (A ^ ( - Ti;) - (ACT - >)2 )' 

where the averages are taken over the 64 points within 
each separate packing. In Fig. [10] we plot S(£cg) an d 
obtain that the inhomogencitics decay as 1/£cg- The 
force fluctuations thus lack an intrinsic length scale. We 
can however collapse the data for S for different pressures 
when they are plotted as function of £cgAz. Therefore, 
the length-scale £* ~ 1/Az, which was observed in the 
inflation response in the previous subsection, can also 
be used to set the coarse graining length one needs to 
obtain a stress field with a particular inhomogeneity S. 
It should be noted that Fig. [TU] also implies that, for a 
fixed coarse graining length, the ensemble size one would 
need to obtain a desired smoothness of response grows as 
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(Az)~ 2 . 



IV. ENERGY MINIMIZATION AND LOCAL 
SLIDING 

In this section we will explore various connections be- 
tween the local displacement field, global energy mini- 
mization, and the scaling of the elastic moduli. In sec- 
tion IIV Al we present simple scaling arguments for the 
typical magnitude of u±, uu and their ratio as function 
of the distance to the jamming point. In section HV Bl we 
employ the displacement angle distribution P(a) Q to 
verify the latter scaling prediction, and find that both un- 
der shear and under compression, the response of jammed 
packings becomes increasingly non-affine near the jam- 
ming point. In section PV CI we test the scaling predic- 
tions for the probability distributions P(uu) and P(wj_), 
and find that u± diverges near the jamming point, as 
predicted. However, while our simple scaling arguments 
capture the behavior of systems under shear in detail, 
they do not quantitatively capture the scaling of -P(it||) 
and P(u±) for jammed packings under compression. 



A. Simple scaling arguments 

We now present a set of simple scaling relations con- 
necting the local displacement field to the scaling of 
the elastic moduli and global energy minimization. The 
starting point for our discussion is the energy expansion, 
equation ([4]), which is a function of the relative displace- 
ments of particles in contact, and which for our power 
law interparticle potential reads (56j: 

^ E = iE^fei-f • ( 22 ) 

(id) V J 

As explained in section ITBl fey is the stiffness d 2 Vij /dr 2 ^, 
where Vy denotes the interparticle potential, Sij is the 
dimensionless overlap 1 — ry/(i£j + Rj), and tin (u±_) 
denotes the parallel (perpendicular) component of the 
relative displacement Uy, with respect to the vector re- 
connecting the centers of the contacting particles. 

For the response to a global shear or compression, AE 
is also related to the elastic moduli, through 

AE ~ Ce 2 ~ i ^2 % («jU - |V' 2 .,,) > (23) 

where C refers to G and K for shear and compressive 
deformations, respectively, and e is the applied strain. 
Note that the first term between brackets is strictly pos- 
itive and the second term is strictly negative. Because 
the packings were constructed to be at an energy mini- 
mum, this implies that the second term cannot dominate 
over the first one, and hence the scaling of typical values 



of tiii is connected to that of the corresponding modulus. 
Inserting K/k ~ 6° and G/k ~ 6 1/2 0, we obtain 

it|| - e<S 1/4 for shear (24) 
uii ~ eS a for compression, (25) 

where symbols without zj-indices refer to typical or av- 
erage values of the respective quantities. Note that by 
dividing out the stiffness k we obtained exponents that 
are valid for both harmonic and hertzian interactions, 
and in fact should be valid for all potentials of the form 

V ~ S a - . 

n 

The expected scaling for u±, the amount by which 
particles in contact slide past each other, is more sub- 
tle. Because of the negative prefactor in Eq. (|22"|) , energy 
minimization should maximize these u±. Again, they are 
bounded by the fact that these packings are stable: unh- 
and u±,ij arc not the independent variables of the prob- 
lem, as they are coupled through the packing geometry, 
and for stable packings the question really is how close 
to the typical u| the typical 8u 2 ± can get. The best the 
system could do in order to minimize the change in en- 
ergy is to make them of the same order, which suggests 
that 

^l~V6. (26) 
u± 

In particular, this means that the system is always close 
to a buckling instability, so that compressing it would in- 
evitably lead to the formation of more contacts to resta- 
bilize the system Q. For the displacements in low en- 
ergy eigenmodes, Eq. ([26]) has recently been derived [43| 
- more intuitively (|26|) can be understood as arising 
from the elastic distortion of floppy modes on a scale 
I* ~ 1/ Az ~ 1/ \f5 [3l[ . Generically, one may expect this 
property to carry over to the displacements in response 
to external perturbations, but this is by no means guar- 
anteed. Indeed, we shall find that Eq. |26|) and the corre- 
sponding predictions for the scaling of u± break down in 
case of a global compression, while they work very well 
for the shear response. In a related paper, we studied 
this issue in depth in the context of networks of harmonic 
springs [3lj . 

It is important to understand the relation between 
Eq. (f2"2"]l and the excess coordination number Az = 
z — Zi SQ . First, note that for z < Zj SO , deformation fields 
exist for which AE — — the so called floppy modes 
(see Appendix [U]). Generically, right at the jamming 
transition (z = z; so ), there are just enough degrees of 
freedom to set all individual tin y equal to zero for the 
type of response considered, in other words, to have side- 
ways sliding motion for all of the bonds at jamming. For 
hyperstatic packings with z > Zj so , the My cannot all be 
set to zero anymore, because there are more of them than 
the number of degrees of freedom dN. In addition, the N 
rf-dimensional coordinates are sufficiently constrained so 
that displacements which make the negative terms in AE 
outbalance the positive ones are forbidden, and AE > 0. 
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A thorough analysis of ([26)) involves the scaling of the 
distribution of the ratios uh ij/u± t ij, since the average of 
uj_ and uu over all contacts may be dominated by large 
and rare fluctuations. To quantify the relative amount 
of sliding and deformation, we therefore introduce the 
displacement angle otij , defined as the angle between u,j 
and Vij, or, 

tana,, = ^ . (27) 

U\\,ij 

In section llV Bl below we will present the probability den- 
sities P(a) for shear and compressive deformations as a 
function of the distance to point J, and test the predic- 
tion given by Eq. ([2"6"| in detail. 

In section IIV CI we will study the scaling of the dis- 
tributions of Mil to test the predictions of Eqs. (|24I25[) . 
We will also present the the distributions of u± under 
shear and compression, and explore to what extent these 
saturate the stability bound given by uj_ ~ uu /y/6. As 
alluded to before, we will show that for the response to 
shear, u± follows the prediction from combining Eqs. ([24]) 
and (|26[) . namely 

u± ~ e<r 1/4 for shear. (28) 

The corresponding scaling for the response to compres- 
sion would be, from Eqs. (|25[) and (|2G|) . 

u± rsj eS^ 1 / 2 for compression, (29) 

but we will show that the numerical results do not sup- 
port this — in the case of compression the two terms in 
the expansion of AE do not become of the same order. 

Note that Eq. ([2"6"| predicts that the response of the 
system becomes strongly non-affine near point J, and 
that Eqs. (|28I29[) predict that typical distances by which 
particles slide past each other diverge as we approach the 
jamming transition. In finite systems, this divergence of 
the local displacements is limited by finite size effects. We 
discuss the corresponding finite size scaling and crossover 
in Appendix [D] 

B. The displacement field and displacement angle 
distribution 

We extract the displacement fields from our linear re- 
sponse calculations of the periodic packings of 1000 par- 
ticles. Examples of these, for both compression and shear 
and for three different pressures, are shown in Fig. Ilia . c. 
and the corresponding displacement angles a are shown 
in Fig.HIt,d. 

For high pressure, far from point J, the displacement 
fields are quite smooth and close to affine deformations, 
i.e., u = e(yx + xy) for shear and u = — e(xx + yy) for 
compression. When point J is approached, the displace- 
ment fields clearly become more disordered. The dis- 
placements become increasingly non-affine, and organize 



in vortex-like structures. Similar structures have been 
observed in various experimental and numerical stud- 
ies of disordered systems close to the jamming transi- 
tion [n, m b m m, S3, si . 

In principle, this non-affine behavior can be analyzed 
by subtracting the appropriate affine displacement field 
from the observed displacement field [H, S|| S3 • By do- 
ing so, however, valuable information on the correlation 
between and u»j is lost. Therefore, we analyze the 
nature of the non-affinity in terms of the displacement 
angles (a), where a is defined via tan c^j = ^p^- (see 

Eq. {27} ). Under an affine compression, all particles get 
closer together along their contact line and P(a) would 
be a (5-peak at a = ir. For a purely affine shear dis- 
placement field in an isotropic system, P(a) is flat [57] ]. 
Finally, close to point J where <5 — > 0, Eq. ([26|) predicts 
that Ux,ij /u\\ij diverges, i.e., that P(a) approaches a 5- 
peak at a = ir/2 (see also Appendix ID]) . 

The values of a in Figs. [TTbl and [TTHl are what is 
expected for an affine displacement. One sees that close 
to jamming, more and more bonds have a close to n/2 
(corresponding to the bright yellow dots in Figs. lllb 2 and 
ITTH2 and in particular Figs. [TTb3 and ITTH3) . indicating 
that for those displacement fields uu <C as predicted 
in Eq. f26]) . 

For the rest of the paper we will ignore the spatial 
organization of the displacement angles a, and focus on 
their probability densities P(a). In Fig. [12] we show these 
probability distributions for comprcssional and shear de- 
formations for a range of pressures. For large pressures, 
P(a) reflects the affine limits discussed above, while upon 
lowering the pressure, P(a) evolves to exhibit a strong 
peak at a = ir/2, as predicted by Eq. P5[) . In the case of 
the shear deformation, we find that the peak in P(a) can 
be well approximated by a Lorentzian, and we fit P(a) 

to [E3 

P(a) ~ ■ (30) 

n w i 4. ( a -*l 2 \ 

Note that the width if should vary as Un/u±_ close to the 
jamming transition, because \otij — 7r/2] ~ MiUj/itj^jj if 

u \\,ij ^ U -L,»j 

In good approximation, the resulting widths satisfy the 
scaling relation 

w ~ <5°- 55±0 - 05 , (31) 

as is shown in the inset of Fig. fT2a . This scaling has 
been used to construct the colored surface in Fig. [T2"a 
from Eq. (|30|) . and is consistent with the prediction of 
Eq. (EH), that w ~ VS. 

For the case of compression the situation is more com- 
plex as can be seen from Fig. [T2"b . Since the distribution 
here does not appear to be governed by a single scale 
and we do not have a fitting function, we have used the 
full width at half maximum of P(a) as w. The value of 
this width is shown in the inset of Fig. \T2b . It scales as 
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FIG. 11: (color online) Examples of displacement fields in linear response to shear (a) and compression (c) for three pressures, 
p = 3 • 10 _2 ,5 ■ 10 _4 ,5 ■ 10 -6 , from top to bottom. In (b) and (d) the corresponding displacement angles are depicted, where 
bonds are marked with a dot, on a color scale which runs from red (a = 0), via yellow (a — n/2) to green (a = n) — see 
Eq. (|27p for the definition of a. Clearly, the highest pressure packing displays almost affine displacements, while at lower 
pressures the response becomes increasingly non-afnne. 



w ~ (5 - 44 ! 5 ], which is reasonably consistent with Eq. (f!?fj)) . 
However, the presence of a clear shoulder for a ~ 7r is a 
strong indication a simple one-parameter scaling actually 
does not hold (3l| . 

The shear data shown in Fig. [T^] are thus in agreement 
with the prediction from the simple scaling argument 
that near point J, Un/ux ~ \/5 (Eq. ([26]) ). Hence, we 
conclude that the displacements in response to shear do 
lead to an approximate balance between the two terms 
in AE (Eq. (p2|) ). The closer one gets to point J, the 
more strongly non-afnne the deformation field becomes, 
in other words, the more the mechanical response of the 
system is different from that of a homogeneous, isotropic 
elastic material. For compressive deformations, the sit- 
uation is somewhat more complicated: the development 
of a peak at a w it/ 2 still signifies the same tendency 
to non-amne sideways sliding motion with an amplitude 



that might be described by a balance of the two terms 
in the energy expression, but the remainder of a signifi- 
cant shoulder at large angles a hints at the fact that the 
response can not be understood completely in terms of 
this balance or a simple one-parameter scaling. For more 
details, see Ref. (3l| . 



C. Scaling of local displacements 

As we discussed in the introduction to this section, by 
linking the scalings of the local relative displacements to 
the elastic moduli, we can predict scaling relations for u\\ , 
and, when the two terms in the energy expansion are of 
the same order, for u±, both for compressive and shear 
deformations. We now test these predictions by plotting 
the distributions of the rescaled displacements P(u\\d n ^) 
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FIG. 12: (color online) The displacement angle distribution 
for various pressures, obtained by averaging over 100 packings 
for each pressure, (a) Under shear P(a) evolves from nearly 
flat (highest pressure) to very sharply peaked at a = n/2 
(lowest pressure). The solid lines represent numerical data, 
the gridded surface Eq. (|30|l . Inset: The scaling of the width 
of the peak as a function of typical overlap 5. The dashed 
line indicates w ~ <5 0,55 . (b) Under compression P(a) has a 
peak at it for high pressures and develops a peak at tt/2 at 
low pressures. Inset: The scaling of the width of the peak 
as a function of typical overlap S. The dashed line indicates 
w~S 0A4 . 



and P(u±5 n± ), where the power law indices n follow from 
Eqs. (|24I25I28I29|) . and check for data collapse. 

The appropriately rescaled distributions of u\\ and hi 
are shown in Fig. [13] for the case of shear deformations. 
Note that we normalize the displacements by the strain 
e [59|. As is clear from Fig.[T31 the scaling collapse is very 
good for both components of the displacement field u^ , 
fully confirming the simple prediction for uii presented 
in Eq. and the prediction for u± in Eq. that 

follows from the assumption that the two terms in the 
energy expansion are of the same order. 

For the displacements under compression we show scal- 
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FIG. 13: Collapse of the probability densities of u± (a) and 
My (b) for a shear deformation. The axes have been rescaled 
according to the prediction of Eq. (|24p and (|28[). Higher pres- 
sures have been omitted from the analysis. 
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FIG. 14: Rescaled probability distributions of tii and u\\ for 
compressional deformations, (a-b) Rescaling of P(u±) and 
P(m||) with the exponents predicted by Eqs. (|25|1 and (|29p 
produces a poor collapse, indicative of the special nature of 
compression of jammed packings, (c-d) When the scaling ex- 
ponents are adjusted by hand, from 0.5 to 0.3 (c), and from 
to 0.05 (d), a reasonable collapse can be obtained. As before, 
in all panels the distribitions for the highest pressure have 
been omitted. 



ing collapses in l ':u. ITTI (ill . In Fig.[T4"k-b we show the dis- 
tributions P(u_l<5 - 5 ) and P{u\\S°), which Eqs. J2SJ and 
(|29|) predict to collapse. Clearly, the collapse is far worse 
than in the case of jammed packings under shear. Sim- 
ilar to P(a), the distributions P(u±) and P(u\\) do not 
appear to exhibit simple scaling for jammed packings re- 
sponding to compression. The best scaling collapse is 
obtained by adjusting the scaling exponents n by hand, 
but even then, Fig. [T4b-d show, that this adjusted scal- 
ing collapse is still less convincing than for systems under 
shear. The exponent n± = 0.30 is significantly different 
form the predicted value of 0.5, reflecting that the dis- 
placements in response to compression satisfy neither a 
balance of terms in the energy expansion nor pure one- 
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parameter scaling [3 11 ] - 

In conclusion, the displacement angle distribution 
shows the development of a peak in P(a) at a = tt/2, 
both for shear and compressive deformations. The width 
of the peak shrinks as y/S, consistent with the energy 
balance argument that predicts that uu/u± \/S. This 
peak signals that near point J, most particles in contact 
mainly slide past each other, in a manner which helps to 
minimize the changes in elastic energy. However, in the 
case of compression P(a) retains a significant shoulder 
close to a = tt, even close to the jamming transition. 

Both for shear and compressive deformations, the slid- 
ing motion u± diverges near jamming. This suggests that 
corrections due to finite size play a role — for details see 
Appendix [Dj For shear deformations, simple arguments 
predict the scaling of P(u\\) and P(u±) quantitavily cor- 
rect, but for compressive deformations our simple argu- 
ments break down. We believe that this fact can ulti- 
mately be traced back to the special geometry that the 
contact network of a packing has because it is made of 
particles that interact purely repulsively. We have re- 
cently studied this question in the context of harmonic 
networks — for more information sec Rcf . 1311 . 



SUMMARY AND OUTLOOK 



We have shown, by means of linear response calcula- 
tions, how the applicability of elasticity theory to disor- 
dered packings of frictionless spheres is related to the dis- 
tance to the jamming transition. Averaging the response 
to a local perturbation over an ensemble of packings, we 
fitted the stress response to continuum elasticity. The 
resulting clastic moduli arc the same as those obtained 
by calculating the response to a global shear or compres- 
sion, proving the consistent applicability of continuum 
elasticity in an ensemble averaged sense. In single pack- 
ings, we identified a length scale t* ~ l/(Az) up to which 
the response is dominated by local disorder. Since I* di- 
verges as the jamming transition is approached, contin- 
uum elasticity breaks down completely as a description 
of the linear response of individual packings near jam- 
ming. The length scale I* corresponds to the length used 
by Wyart et al. to estimate the number of soft modes in 
systems near jamming |(|. Here we have shown that it 
also represents the coarse graining length needed to ob- 
tain a smooth stress response tensor in a single globally 
deformed packing. 

The displacement fields near jamming arc highly non- 
affine. Wc analyzed these displacements at the grain scale 
through the statistics of local relative displacements of 
neighboring particles. These relative displacements, de- 
composed into U|| along the line of contact, and u± per- 
pendicular to it, govern the elastic behavior of the system 
because according to Eq. (|U) they enter into the change 



m energy as 



AE 



*5> 

(hi) 



We have introduced the displacement angle distribution 
P{a) to easily characterize the non-affinities, via the an- 
gle between the relative displacement vector and the vec- 
tor connecting two neighboring particles. Close to point 
J, these distributions develop a peak near a = n/2, in- 
dicating that the non-affinity is such that particles pre- 
dominantly slide past each other. 

Wc also determined the scaling of typical values of u± 
and My separately. The parallel displacements U|| follow 
a scaling that is consistent with the scaling of the elas- 
tic moduli: they are nearly independent of distance to 
jamming for compression and vanish as S 1 ^ 4 for shear. 
The perpendicular displacements u± diverge in response 
to both shear and compression, consistent with the pre- 
dominance of sliding near jamming. 

There is a surprising and fundamental difference be- 
tween compression and shear. For shear, u±_ simply scales 
as so that both terms in the elastic energy are of 

the same order, and the scaling of P(a) is captured by a 
single parameter. For bulk compression, all of this breaks 
down — u± still diverges, but not with an exponent that 
follows from a simple scaling argument, and there is no 
one-parameter scaling for P{a). The issue of what is spe- 
cial about the compression of packings is discussed in a 
separate paper 31 1 . 

There are many extensions of this framework beyond 
the linear response of frictionless discs in a computer 
simulation, of which the extension to three-dimensional 
spheres is probably the simplest. For finite displace- 
ments, the direct connection between displacements and 
energy is lost. However, in any system of approximately 
spherical particles close to the jamming transition, where 
steric hindrance becomes an important factor in the dy- 
namics, one would expect the displacement angle distri- 
bution to develop a maximum near a = 7r/2, correspond- 
ing to predominantly sliding displacements. These lo- 
cal displacement statistics are accessible experimentally, 
even outside the linear response regime, for example in 
steadily sheared (wet) foams or emulsions, and even in 
systems below t he j am ming density, for example in driven 
granular gases [lij HO, H2|- ^ n the latter case, one 
has to define neighbors (e.g. through a Voronoi tessela- 
tion), and one has to choose the time scale t over which 
to measure the relative displacements so that neighbor- 
ing particles have had the time to interact — P{&) will 
depend on r, and preliminary studies suggest that the 
peak in P(a) is maximal for a finite value of r. Note 
that in dense foams or emulsions, in which the particles 
are continuously interacting, one rather expects the peak 
of P(a) to be maximal for t-*0. 

In systems of ellipsoids [H, [H, [l?], EH or frictional 
particles [HI, EH , rotations and torques come into play, 
and the expression for the elastic energy will be more 
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complicated than Eq. ([4J. However, in linear response 
the elastic energy will still be a function of the local rel- 
ative changes of the particle coordinates, and one can 
still check for the occurrence of combinations of local 
displacements that lead to a small or vanishing change 
in energy. For example, in packings of frictional spheres, 
the locally floppy displacement is for contacting particles 
to roll without slipping. For ellipses "zero energy" lo- 
cal displacements can be identified similarly. We expect 
that such low energy displacements should dominate the 
statistics when the system is close to the relevant isostatic 
limit, and we expect that appropriate generalizations of 
P{cx) for such systems may serve as a useful characteri- 
zation of the non-afhnc deformation field. 
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APPENDIX A: PACKING GENERATION 

1. The Hertzian interaction 

The Hertzian interaction is taken from Refs. [321I36T] in 
terms of a force law as 



3/2 



(Al) 



where Rij is the reduced radius of the pair of particles, 
defined by 



1 



1 

Ri 



1 



and E*j is a similar combination of the modified Young 
moduli, 



1 



i 

~E* 



1 

~E* 



The modified Young modulus is determined from the 
Young modulus E and the Poisson ratio v of the ma- 
terial that the grains are made of: 



E* = 



E 



Now let us determine a closed expression for the prefac- 
tor £y of the interaction as we have used it throughout 
the paper. Differentiating the interaction in equation ([1]) 
gives 



f = V-L ( 1 - liL 



3/2 



with dij = Ri 
2| yields 



Rj 



(A2) 

Equating the equations (|A1[) and 



j5/2 



Our unit of pressure is E* (and all grains arc made from 
the same material) , so that in our reduced units E*j = 5 
and 

(A3) 

With the radii of the particles drawn from a flat distri- 
bution between 0.4 < R < 0.6, the value of eij can vary 
between 0.068 and 0.230. To check how large the effect 
of this varying e is, we have performed some simulations 
using a constant e = yg, the value it has for R = 0.5, and 
did not observe significant deviations. This justifies com- 
paring our results to other simulations done at constant 
el. 



APPENDIX B: 2D ELASTICITY AND POINT 
RESPONSE 

The two-dimensional version of linear elasticity has the 
same form for Hooke's law as the well-known 3D ver- 
sion 0, n 

a al3 = 2fiu a p + Xu^5 a /3 , 

where u a p is the strain tensor, S a p is the Kronecker delta, 
and we use the summation convention. A and fi are 
known as the Lame coefficients. All dependence on di- 
mensionality basically arises from the S a p, through the 
fact that the trace of a a fj reads 

a aa = (2/i + d\)u aa . 

The shear modulus G is just the same as the Lame co- 
efficient n, independent of d. The definition of the bulk 
modulus K is the resistance to change of "volume" , 



K 



P 



V/V 



In an isotropic material the stress tensor corresponding 
to a uniform compression is a a p = p5 a p. Its trace is 
therefore o~ aa = pd. The trace of the strain tensor is 
u aa = V/Vq, so that in d = 2, we now have for the 
moduli 



l-I/ 2 



G 
K 



li ■ 

0~n 



du aa 



, 2 



(Bl) 
(B2) 
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The Navier-Cauchy equation, which follows from insert- 
ing Hooke's law and the definition of the strain tensor 
into the equation of mechanical equilibrium d@cr a p = 
— /„ xt , reads (in arbitrary dimension), 

/iAu + (A + //)V(V ■ u) = -f oxt , 

where A denotes the Laplacian operator. In d = 2 the 
coefficients happen to become simply 



GAu + KV(V ■ u) = -f cxt , 
which is equation (fT4| . 



(B3) 



1. Response to a point force 

For a rectangular packing with periodic x-boundaries 
and hard walls on top and bottom, the Navier-Cauchy 
equation (|B3|) can be solved in terms of a Fourier series. 
For a point force in the middle, taken to be the origin of 
the coordinate system, we expand: 



^ A lm sir 



lm 

U V = 

lm 

f x = o 



irmy 



2nlx 

Bi m cos I — — I cos 



(B5) 
(B6) 

/„ =J2F lm co S ( 2 ^)cos(^), (B7) 



L, 



where 



F, 



I in 



L 2 L m odd 
to even 



(B8) 



so f y represents a 8- force in the origin of unit weight. The 
solution is to be compared to granular stress fields which 
have been coarse grained using a Gaussian. Using the 
same Gaussian instead of a <5-function for f y here leads to 
a better fit and to improved convergence of the resulting 
Fourier series for the stress tensor. Since convolution 
in real space is just multiplication in Fourier space this 
amounts to multiplying Fi m by 



L T L. 



■ cxp 



,^2 2 
-7T <T 



2P 
LI 



2LI 



Inserting the expansions into equation (|B3|) we can 
solve for Ai m and Bi m . Solutions for the strain tensor 
and stress tensor are then obtained by taking the appro- 
priate linear combinations of derivatives of u x and u v . 
The resulting stress tensor only depends on the Poisson 
ratio v = [K — G)/{K + G). A separate fitting pro- 
cedure for the displacement field then gives K + G, us- 
ing equation (|16p and the determination of the moduli is 
complete. 



APPENDIX C: FLOPPY MODES 

Suppose we are at point J, so 8 = 0, z = 4, and the 
number of contacts equals 27V, which is equal to the num- 
ber of independent displacement components M, >a . When 
one contact is removed, we can in principle write down a 
displacement field for which all My y = 0. This displace- 
ment field has energy zero: it is a floppy mode [y, [37| . 
This is not merely an artifact of the expansion in u\\ and 
u±: the counting of variables and equations in princi- 
ple allows to write down a displacement field for which 
the unexpanded Ary from Eq. ([3]) are zero and hence 
the change in energy is strictly zero. For a floppy mode 
distortion, AE is identically zero to all orders in the dis- 
tortion, even if the initial configuration did not obey force 
balance. Let us check this up to 0{u\_) in the expansion 
of AE. To see this, we include the term linear in un, 
which we left out of the original equation ^ because it 
vanishes when summed over all contacts: 



AE = 



fij u \\,ij + rfkij 



I, 



u ll,y 



13 u 2 , 



(B4) For a floppy mode, Ary = 0, and this implies that Uiuj 



— Uj_ ij/(2rij) + 0(uj_). Hence, the first and third term 
in the expansion of AE cancel without having to sum 
over all contacts and we are left with AE = 0(u\). This 
even holds if the initial system would not satisfy force 
balance 1611 . 



APPENDIX D: FINITE SIZE EFFECTS 

If the system is sheared or compressed, local relative 
displacements diverge upon approaching the jamming 
point, according to the analysis in section IIV CI How- 
ever, the relative displacement between a given particle 
i and its image in the neighboring unit cell of the peri- 
odic packing is given by the imposed boundary condition 
— hence there is coupling between the magnitude of the 
local relative displacements and the system size. In ad- 
dition, we have extracted the length scale i* from the 
elastic behavior of our system, and should consider what 
happens when this length scale becomes of the order of 
the system size. We tentatively suggest the following pic- 
ture for the implications of this coupling. 

Let us focus on shear deformations because the scal- 
ings are cleanest in that case. The x-component of the 
relative displacement between a particle at (xi,yi) and 
its periodic image at (xj, yi + L) can be written as a sum 
along a path from the particle to its image: 



path 



eL 



(Dl) 



where e is the applied shear strain. If the system is small 
compared to the length scale £*, we do not expect to 
see any elastic behavior, and the displacement field is 
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FIG. 15: (a) Scaling of typical Wj 

p. The dashed line indicates u± ~ p -1 ^ 6 , to which all large 
packings stay close, except for the bottom-right points which 
are very far from point J (p — 3 • 10~ 2 ). The system size de- 
creases from N = 10 000 (diamonds), N = 1000 (triangles), 
N — 300 (squares) to N — 100 (stars). The curve seems 
to level off at small p for the small systems, as predicted 
in equation (|D3|I . (b) The same data, plotted according to 
the finite size scaling prediction of equation (|D7[) . The plus 
signs indicate data taken from the linear response to shear 
of bidisperse packings generated using Conjugate Gradient 
minimization (courtesy N. Xu), which allows to probe much 
lower pressures. The two rightmost points in this graph rep- 
resent series of packings which are essentially isostatic, and 
hence have a very large and inaccurate £* . The dashed line 
represents the predicted behavior for small I* /L. 



dominated by the local packing disorder. In this case, 
we might just as well view the displacements in the in- 
direction as we go along the path as a random walk. This 
walk consists of approximately L steps of size roughly u±, 
so that 



path 



ij,x 



U± 



Vl 



Equating (|Dlj) and (|D2|) yields 

U± ~ eV~L for L <C 



(D2) 



(D3) 



On the other hand, large systems can be treated as a 
continuum down to the scale £*, so if we want to apply 



this model to a large system, we have to consider the 
"random walk" in a subsystem of size £* , and assume the 
globally applied shear deformation scales down affinely to 
that scale. Therefore the boundary condition becomes 



E 



et 



subsystem 

and the random walk in the subsystem gives 



Ui - 

subsystem 



so that finally 



Uj_ 



for L > 



(D4) 



(D5) 



(D6) 



Hence the argument reproduces our result for the scaling 
of u± under shear, equation (f28|) . because £* ~ <5 _1 / 2 . 
Furthermore, it predicts that u± should saturate around 
the value it has when L = £* when approaching the jam- 
ming transition in a finite system. The two limits can be 
connected by a finite size scaling function 



Uj_ 



(D7) 



where f(x) — ► const as x 3> 1 and f(x) ~ x as x <C 1. 
We do not have enough data on systems where L -C £* 
to prove this scaling prediction with firm confidence, but 
the data we have is at least consistent with it. In Fig. [TBI 
we show the typical values of uj_ as a function of p for 
various system sizes, taking £* = 6/Az in accord with 
(f2"Tj) . Figure [T5*b shows the scaling collapse correspond- 
ing to equation (|D7[) . The two rightmost points in this 
plot have z m 4, so £* is very large and very inaccurate, 
but nevertheless some levelling off of the curve can be ob- 
served. More data using packings between 100 and 1000 
particles at 4.005 < z < 4.02 could shed more light on 
the scaling behavior. 
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The point which deviates from the scaling of G is ob- 
tained for the lowest pressure from the point response. 
We believe this large deviation to be due to the fitting 
procedure: we fit the Poisson ratio v from the stress field 
Aa a p(r). As p — > 0, we observe v — > 1, its maximum 
value in two dimensions. In this limit, the fitting of the 
stress tensor becomes less accurate for numerical reasons, 
which is made even worse by the fact that the data be- 
come more noisy at low p. In the bulk response we see a 
similar effect: the fluctuations around the average values 
are much larger at lower pressure. 
See page 21 of Ref. 

For a general potential Vij ~ 5"j the factor 2/3 in (|22|) 
becomes 1/(q — 1). 

This is easily seen in the case of a simple shear in the 
x-direction: All Ui are in the z-direction, hence so are 
the Uij, while the ry are isotropic. 

Because of the 7r-periodic nature of the angle a, we add 
copies of the peak at —n/2 and at 3tv/2. This improves 
the fit for larger pressures. 

In addition, it is natural to eliminate systematic depen- 
dence of My and ui on the particle radii, and we have di- 
vided each u± t ij and Miuj by the actual length nj of the 
contact — this would also allow to compare the results 
to the case of affine displacements, which have local w±,y 
and Uii y which are proportional to r^, but we have al- 
ready seen that the displacements are strongly non-affine, 
so we do not show the corresponding plot here. 
Note that because the distribution of u\\ under compres- 
sion is not centered around zero, we are collapsing the 
average of the distributions, not the width. 
This holds in general for floppy modes in systems with 
radial interactions between the particles. 



